## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.2 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.3 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
##
## Attachement du package : 'plotly'
##
##
## L'objet suivant est masqué depuis 'package:ggplot2':
##
## last_plot
##
##
## L'objet suivant est masqué depuis 'package:stats':
##
## filter
##
##
## L'objet suivant est masqué depuis 'package:graphics':
##
## layout
main_theme = theme_bw()+
theme(line = element_blank(),
axis.line = element_line(colour = "black"),
panel.border = element_blank(),
axis.ticks = element_line(colour = "black"),
axis.text.x = element_text(colour = "black", size=22),
axis.text.y = element_text(colour = "black", size=22),
legend.title = element_text(colour = "black", size=20),
legend.title.align=0.5,
legend.text = element_text(colour = "black", size=18),
axis.title=element_text(size=28))
#Parametres
i = 20 # Nombre d'especes
Xmin = 0
Xmax = 2
X0 = (Xmax - Xmin)/2
K0 = 1
lambda = 1
sigma = 1 #etalement de la competition
r = 1
############ Fonctions #############
K_cap = function(x, K0_ = K0, lambda_= lambda, x0 = X0){
(K0_ - lambda_*(x-x0)**2)
}
a = function(x1,x2, sigma_a = sigma){
exp( (-1/2*(x2-x1)**2)/sigma_a**2 )
}
fitness = function(x1, x2, r_ = r ){
r_*(1-a(x1, x2)*(K_cap(x1)/K_cap(x2)))
}
## Set range and domain of plot
x1_interval <- seq(0.1, 1.9, length.out = 25);
x2_interval <- seq(0.1, 1.9, length.out = 25);
## Interpolate surface
z <- outer(x1_interval,x2_interval,
FUN = fitness)
p <- persp(x1_interval,x2_interval,z, theta = 30, phi = 20,
col = "lightblue", shade = 0.4, ticktype = "detailed")
p
## [,1] [,2] [,3] [,4]
## [1,] 9.622504e-01 -0.1900112 0.5220515 -0.5220515
## [2,] 5.555556e-01 0.3291090 -0.9042196 0.9042196
## [3,] -1.567806e-17 0.4812177 0.1751489 -0.1751489
## [4,] -1.517806e+00 0.3651056 -2.1663676 3.1663676
plot_ly() %>% add_surface(x = ~x1_interval, y = ~x2_interval, z = ~z)
## Set range and domain of plot
x1_interval <- seq(0.1, 1.9, length.out = 500);
x2_interval <- seq(0.1, 1.9, length.out = 500);
## Interpolate surface
z <- outer(x1_interval,x2_interval,
FUN = fitness)
couleurs <- ifelse(z > 0, "lightblue", "black")
matrice_couleurs = couleurs
df <- data.frame(x = rep(1:nrow(matrice_couleurs), each = ncol(matrice_couleurs)),
y = rep(1:ncol(matrice_couleurs), times = nrow(matrice_couleurs)),
couleur = as.vector(matrice_couleurs))
ggplot(df, aes(x = x, y = y, fill = couleur)) +
geom_tile() +
scale_fill_identity() +
labs(title = "PIP") +
theme_minimal()+
labs(fill = "Légende\nNoir = Négatif\nBleu = Positif")
K_cap = function (X){
max(0, K0 - lambda*(X - X0)^2) + 10^(-9) # On rajoute 10^-9 pour ne pas avoir de pb de division par 0
}
sigma = 0.3
############### Simuler LV pour un nombre arbitraire d'espece ############
#Conditions initiales
x = seq(from = 0, to = 2, length.out = i) # Valeurs des phenotypes
k = c()
for (i in 1:length(x)){
k = c(c(k),c(K_cap(x[i])))
}
# K de chaque espece
n = rep(1/i,i) #densite de pop initiale
# On met sous forme de matrice
N0 = t(as.matrix(n))
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)
X = matrix(rep(x,i), ncol = i, nrow = i, byrow = TRUE)
#Calcul de la matrice M
D = t(X) - X
A = exp(-0.5*D^2/(sigma^2))
M = A/t(K)
#Resolution du systeme d'ODE
model <- function(t,N,sigma){
dN <- r * N * (1 - N %*% M)
return(list(dN))
}
ode <-ode(N0,0:100,model,parms = sigma)
# Representation
#barplot(ode[100,1:i+1])
ode <- as.data.frame(ode)
#colnames(ode) <- c("t",'x1','x2','x3')
ode_plot = ode%>%
pivot_longer(-time, values_to = "dens", names_to = "sp_ID")
ggplot(ode_plot)+
geom_line(aes(time, dens, col = sp_ID))
ode_plot$sp_ID <- as.numeric(ode_plot$sp_ID)
ode_plot$trait <- x[as.numeric(ode_plot$sp_ID)]
ggplot(ode_plot)+
geom_raster(aes(trait, time, fill=dens))+
scale_fill_gradient2(low = "white" ,
high = "red")
######## Nombre d'espece qui persiste en fonction de sigma (pour l'instant ne marche pas..) ##############
seuil <- 0.05
nb_especes <- c()
for (sigma in seq(0.3,5,by = 0.1)){
ode <-ode(N0,0:100,model,parms = sigma)
nb = 0
for (m in tail(ode,1)[,1:i+1]){
if (m>seuil){
nb = nb + 1
}
}
nb_especes <- c(c(nb_especes),c(nb))
}
plot(seq(0.3,5,by = 0.1),nb_especes)
######## Implementation de mutation ##########
#parametres
p = 0.1 #taux de la pop qui mute
T = 10000
t = 200
n = rep(0,i)
n[i/6] <- 1 # densites de pop initiale (une seule espece)
N0 <- t(as.matrix(n))
#Initialisation des valeures de K
k = c()
for (i in 1:length(x)){
k = c(c(k),c(K_cap(x[i])))
}
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)
M = A/t(K)
ode <-ode(N0,1:t,model,parms = sigma)
ode_mutation <- ode
#trouver un moyen de faire muter tout le monde... (peut etre metre un seuil de densité min..)
for (z in 2:(T/t)){
N0 <- t(as.matrix(tail(ode,1)[,1:i+1]))
a = runif(1)
index = which.max(N0)
if (a<0.5){ #mutation vers la gauche
N0[max(0,index-1)] <- p*N0[index]
N0[index] <- (1-p)*N0[index]
} else { #mutation vers la droite
N0[min(i,index+1)] <- p*N0[index]
N0[index] <- (1-p)*N0[index]
}
ode <-ode(N0,(z*t):(z*t+t),model,parms = sigma)
ode_mutation <- rbind(ode_mutation,ode)
}
ode_mutation_plot = as.data.frame(ode_mutation)%>%
pivot_longer(-time, values_to = "dens", names_to = "sp_ID")
ode_mutation_plot$sp_ID <- as.numeric(ode_mutation_plot$sp_ID)
ode_mutation_plot$trait <- x[as.numeric(ode_mutation_plot$sp_ID)]
ggplot(ode_mutation_plot)+
geom_raster(aes(trait, time, fill=dens))+
scale_fill_gradient2(low = "white" ,
high = "red") +
main_theme
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
# ggplot(ode_mutation_plot)+
# geom_tile(aes(trait, time, fill=dens))+
# scale_fill_gradient2(low = "white" ,
# high = "red") +
# main_theme